Local adaptation to climate inferred from intraspecific variation in plant functional traits along a latitudinal gradient

Results of this study show how the morphology and physiology of an endemic Australian plant respond to local climate and soil profiles across a ~160 km latitudinal gradient. Climatic factors accounted for more of the variation in functional traits than soil, suggesting that trait shifts might occur with climate change.


Introduction
Anthropogenic activities are causing widespread and rapid changes to the earth's climate, with weather extremes, such as frequent heatwaves, prolonged severe droughts and variability of precipitation events increasing (IPCC, 2021).These climate events are expected to have major consequences for the fitness and persistence of numerous species, with risks of severe spatial displacements and extinctions (Fletcher et al., 2020;Miller et al., 2020).Impacts associated with climate change are some of the most challenging threats to manage, given the global scale and the uncertainties of the responses of species and ecosystems to multifaceted drivers of environmental change (Millar et al., 2007;Bellard et al., 2012;Dube et al., 2018).Additionally, there is growing need to understand how populations vary in their acclimation and adaptation to environmental change, and the risk of local populations not possessing the plasticity of phenotypic traits to effectively cope with changing conditions, or the genetic variability to adapt to them (Breed et al., 2018;Harrison, 2021).Examining how populations vary in their responses to interacting environmental factors delivers key insights into the drivers of plant form and function, as well as the patterns and processes potentially underpinning species resilience or vulnerability to those environmental factors.
Trait-based approaches provide an empirical foundation for explaining patterns of variation and forecasting consequences of environmental change to ecosystems (Bjorkman et al., 2018;Heilmeier, 2019).Functional traits can be considered as any morphological, physiological or phenological attribute of an organism, which can ultimately impact a measure of performance or fitness (Lavorel et al., 2007;Violle et al., 2007).These attributes can include indicators of plant form or function, such as plant size and photosynthetic rate (Cornelissen et al., 2003;Perez-Harguindeguy et al., 2016).Plant functional ecology aims to establish generalisable predictions across biological, temporal and spatial scales (Adler et al., 2014), but calls have been made to better understand relationships between functional traits and the environment with explanatory and predictive power (Shipley et al., 2016;Des Roches et al., 2018).Established trait-based approaches have historically focussed largely on trait variation between species (McGill et al., 2006).However, one of the primary challenges in trait-based ecology is accounting for and identifying the factors driving intraspecific variation in functional traits (Violle et al., 2012).This is grounded in the knowledge that functional traits vary at the individual level and can have cascading effects across biological scales, which ultimately drives the functioning and persistence of populations and communities (Bolnick et al., 2011).
Intraspecific trait variation (ITV) emerges from ontogenetic, phenotypic and genetic differentiation (Leimu and Fischer, 2008;Mason et al., 2013;Henn et al., 2018).The ability to produce different morphological and physiological responses to mitigate against environmental stress can increase the capacity for populations to persist into the future (Nicotra et al., 2010;Araújo et al., 2021).Greater ITV is likely to allow species to cope better with environmental stress through phenotypic plasticity in the short term, potentially promoting genotypic adaptation to future environmental norms in the long term (Nicotra et al., 2010).However, research across several biological scales suggests that ITV is often constrained under environmental stress, such as warming, drying or nutrient limitation (Valladares et al., 2007;Grassein et al., 2010;Gugger et al., 2015;Stotz et al., 2021).Reduced intraspecific variation may also be a result of reduced genetic diversity, environmental filtering and strong stabilizing selection for certain traits, potentially leading to individuals being highly adapted to their local environment and improving their performance within a particular suite of conditions (Bradshaw, 1965;Schlichting, 1986;Schneider, 2022).Therefore, measuring the environmental factors that may drive trait variation can provide insight into how individuals are responding to current conditions, and how they may perform in the face of environmental change (Arnold et al., 2019).
Selection pressures along environmental gradients often give rise to morphophysiological traits that show marked and predictable responses that are functionally linked to climatic and edaphic conditions along the gradient (Reich et al., 2003;Dwyer and Laughlin, 2017;Kandlikar et al., 2022).Climate factors, such as temperature and precipitation, impact the distribution and persistence of species through effects on cellular processes, reproduction and survival (Suzuki et al., 2014).Similarly, edaphic factors, such as nutrient availability, acidity, salinity and metal toxicity, can also present critical limitations to plant physiology and growth, via impaired metabolic processes, enzyme activity and membrane stability (Ahanger et al., 2016).Indeed, climatic and edaphic factors can also interact.For example, shifts in temperature and precipitation regimes can drive microbial activity and, by extension, decomposition rates, soil nutrient turnover, organic carbon (OC) stocks and nutrient availability (Lloyd and Taylor, 1994;Figueiredo et al., 2018;Jansson and Hofmockel, 2020).More specifically, warming and drying can suppress microbial activity (Allison and Treseder, 2008), with nutrient limitations often linked to drier, lower rainfall systems (Moreno-Jiménez et al., 2019).Despite their pervasive independent and interactive effects, there are few studies that simultaneously consider both climate and soil on intraspecific variation in plant functional traits in situ (Maire et al., 2015;Simpson et al., 2016;Joswig et al., 2022).Nevertheless, understanding plant responses to local environmental conditions is crucial for guiding the conservation of extant populations, as well as the restoration of those that are degraded, which can maintain species functional diversity and ultimately biodiversity in the face of rapid environmental change.
Increased temperatures and declining rainfall in southwest Australia risk the 'transition or collapse' of the vulnerable Northern Jarrah Forest (NJF) ecosystems (Lawrence et al., 2022).Therefore, this study assessed the physiological and morphological variation within and between five spatially distinct populations of the NJF endemic perennial, Stylidium hispidum (Lindl.),across a ∼160-km transect that spanned the latitudinal extent of the species distribution.We examined population differentiation in trait means and magnitude of trait variability along the latitudinal gradient and determined what proportion of morphophysiological variation was accounted for by the independent and combined effects of climate and soils.We set out to test whether northern populations would display reduced trait variability due to the hotter and drier conditions increasing the likelihood of abiotic stress (Fig. 1), as ITV is often constrained under environmental stress (Valladares et al., 2007;Grassein et al., 2010;Gugger et al., 2015;Stotz et al., 2021).Moreover, we expected trait means to vary along the gradient (Fig. 1), such that northern populations, exposed to warmer and drier conditions, would reflect water conservation strategies such as reduced stomatal conductance, greater water use efficiency (WUE), smaller plant sizes and fewer flowers (Descamps et al., 2021;Valliere et al., 2021), in contrast to southern populations.Given the southern populations are exposed to higher rainfall and cooler climate conditions, we expected that they experience higher nutrients availabilities and reflect traits associated with higher productivity (i.e. higher gas exchange rates, larger plant size, increased number of flowers).

Study species
Stylidium (Stylidiaceae), commonly known as trigger plants, is a large genus, comprising over 300 species distributed over a range of ecosystems and climate niches, mostly in Australia, with over 130 species occurring in the Jarrah Forest region of Western Australia (Western Australian Herbarium, 1988).
The flowers of Stylidium are unique and feature a motile column composed of fused staminate and pistillate tissues that can be 'triggered' by touch.The swift movement is a pollination mechanism that deposits or retrieves pollen from insect pollinators such as bee flies (Diptera: Bombyliidae), hoverflies (Diptera: Syrphidae) and bees (Hymenoptera: Apoidea; Armbruster et al., 1994).
Stylidium hispidum Lindl., or the White Butterfly Triggerplant, is a small perennial herb, ∼12-28 cm tall, endemic to the NJF of southwest Western Australia and represents the focal species for this study.The species is a member of the Stylidium piliferum R.Br.complex and is characterized by a basal rosette with one or more floral racemes that emerge and flower during austral spring.Young plants are composed of a single rosette from an unbranched stem with green, linear, incurved leaves that may be tinged with pale maroon tones.Older plants typically feature up to eight tightly clustered rosettes, which are raised above the ground on stilt roots (Lowrie and Kenneally, 2017).The distribution of S. hispidum is restricted to a narrow transect parallel to the Darling Plateau with plants growing in highly weathered, lateritic soils.The genus has mycorrhizal associations (Brundrett and Abbott, 1991); however, it is unknown whether these mycorrhizas play a significant role in nutrient acquisition (Lambers et al., 2018).Genetic analysis of spatially distinct populations of S. hispidum found significant differentiation among populations separated by at least 13-23 km along the latitudinal transect of distribution with 23% of the molecular variance partitioned among populations (Hufford et al., 2016).In addition to the intrinsic conservation value of understanding how trait diversity is distributed across the distribution of a local endemic species, S. hispidum serves as a tractable representative from which meaningful conservation insights may be derived, especially for the many Stylidium that are already listed as 'priority' and 'threatened' species under the Biodiversity Conservation Act 2016 (Department of Biodiversity Conservation and Attractions, 2016).

Study area
Five study sites were selected along a ∼160-km north-south gradient of the NJF along the Darling Range spanning Julimar (−31.484691, 116.166501) in the north to Hoffman (−32.958298, 115.955832) in the south (Fig. 2).Sites were selected based on populations included in genetic analysis conducted by Hufford et al. (2016) (Gardner and Bell, 2007;Koch and Samsa, 2007).Broadly, the region is characterized by a dry Mediterranean climate that presents hot dry summers and cool wet winters.Sites in the north experience a mean annual rainfall of ∼700 mm compared to ∼1200 mm in the south (Fig. 2), while the mean maximum temperatures are approximately 1.5  (Nemchin and Pidgeon, 1997;Anand and Paine, 2002) and comprise well drained, shallow to moderately deep gravelly sands overlaying lateritic duricrust.The nutrient concentrations in this ancient and deeply weathered landscape are typically very low, especially for phosphorus, which limits plant growth despite several adaptations that maximize the efficiencies of phosphorus uptake and use (Lambers et al., 2011;Lambers et al., 2015;Viscarra Rossel and Bui, 2016).
As much as practically possible, sites were selected to feature similar topographies and community compositions.

Climate and edaphic characterization
Climate conditions were quantified for soil microclimate, as well as local atmospheric climate conditions for each site.Hydrological (soil water content at 10-cm depth) and thermal conditions were measured through three TMS moisture and temperature loggers (TMS4, TOMST Ltd, Prague, Czech Republic) positioned in random locations at each site and active for the duration of the study, over October, November and December of 2021.Local climate data (e.g.cumulative rainfall) for each site were also obtained from weather station data available through the Bureau of Meteorology.We extracted mean annual temperature (BIO1, MAT) and mean annual precipitation (BIO12, MAP) data for each site (Fig. 2) from WorldClim version 2.1 at a spatial resolution of 30 arcseconds (Fick and Hijmans, 2017).Atmospheric vapour pressure deficit (VPD a , kPa) and relative humidity (RH%) were captured using a LI-6400XT portable gas exchange system (LI-COR Biosciences, Lincoln, NE, USA) and averaged across months for each site, by taking the 'reference' measures of the atmospheric air that was taken in through the gas exchange system with desiccant column on 'full bypass'.

Plant morphological and physiological traits
Field measurements were conducted in October, November and December 2021.Each month, morphological and physiological traits were measured across at least 10 reproductively mature plants at each site that had grown to a rosette diameter of at least 2 cm and separated from other sampled plants by at least 2 m.Morphological traits measured include rosette size (product of the perpendicular diameters; length, width and height; cm 3 ), number of flowers and plant height (cm).All morphological traits were paired with the physiological measurements from the same plant, and no plant was repeatedly measured.
Physiological performance was characterized by photosynthetic CO 2 assimilation rate (A), stomatal conductance (g s ), transpiration (E) and leaf temperature using a LI-6400XT portable gas exchange system (LI-COR Biosciences, Lincoln, NE, USA) that was equipped with a 6400-15 Extended Reach 1-cm Chamber receiving ambient light and temperature conditions.The size and morphology of Stylidium rosettes limit the use of other larger chambers.Measurements were taken with chamber conditions set to carbon dioxide (CO 2 ) concentrations of 400 μmol CO 2 mol −1 ; RH and temperature were maintained at ambient conditions to reflect site seasonal conditions for climate.Rosette leaves were displayed as fans upon insertion into the leaf chamber, ensuring minimal overlap between leaves.Measurements were taken when stability was reached, determined via visual analysis of real-time response curves accompanied by system stability indicators.At least three measurements were taken per plant.All gas exchange measurements were adjusted by the leaf area, by capturing a bird's eye photograph of the leaves within the chamber and through pixel analysis using ImageJ.After gas exchange measurements were obtained, intrinsic WUE was calculated as A/g s .
Chlorophyll content was measured using a CCM-300 Chlorophyll Content Meter (Opti-Sciences, Hudson, NH, USA) by placing the sensor head directly on the leaf surface of each plant for three independent leaves per plant.Chlorophyll fluorescence (F v /F m ) was captured using a Pocket PEA chlorophyll fluorometer (Hansatech Instruments Ltd, UK) with photosynthetically active radiation (PAR) set to 3500 μmol m −2 s −1 following at least 10 minutes of dark adaptation on leaves.All physiological measurements were constrained to a sampling period between 9:00 am and 12:00 pm, representing the time when the plants were most photosynthetically active, and conducted within the same

Patterns of latitudinal variation
All statistical analyses and ordinations were performed in the R statistical environment (version 4.1.3)using RStudio Version 2022.02.0 (R Development Core Team, 2022).Summary statistics are presented as means ± standard error for site-level differences.
Differences in morphological and physiological traits within (temporal) and between (spatial) sites were assessed using two-way analysis of variance (ANOVA) and post hoc Tukey tests using the anova and TukeyHSD functions in the 'stats' package, respectively.Prior to analysis, data were assessed for homogeneity of variance using the leveneTest function in the 'car' package, and data were log transformed where appropriate to meet the assumptions for the ANOVA.Phenotypic variation was estimated using the quartile coefficient of dispersion (QCD), calculated as [Q3 − Q1] / [Q1 + Q3] for each trait per population for each sampling month.This statistic is suggested to be a robust alternative to the coefficient of variation, which has been found to underestimate intraspecific trait variability, or plasticity, in approximately 50% of cases (Rosas et al., 2019;Hurtado et al., 2020;Yang et al., 2020).
A Bray-Curtis dissimilarity matrix was used for a nonparametric permutational multivariate ANOVA (Anderson, 2005) to identify significant differences between trait assemblages along the latitudinal gradient using the adonis function in the 'vegan' statistical package.Post hoc pairwise comparisons between trait assemblages were conducted to identify the multivariate differences between sites using the pairwise.perm.manovafunction in the 'RVAideMemoire' package (Anderson, 2014;Hervé and Hervé, 2020).
Principal components analysis (PCA) was used to identify the major dimensions of variation in the trait, climate and soil datasets and visualize multidimensional variation among populations and sites using the PCA function in the 'FactoMineR' package (Lê et al., 2008).The PCA was carried out on log-transformed data, and trait correlations with the first two principal components were assessed using the dimdesc function.

Climatic and edaphic associations
Redundancy analysis (RDA) was used to construct an ordination of the variation in trait composition as constrained by the edaphic and climatic variables and identify dimensions of trait variability that correlate with environmental variables using the rda function in the 'vegan' statistical package (Oksanen et al., 2007).Beginning with the intercept-only model, we applied a stepwise selection procedure using the ordistep function to determine the most parsimonious subset of explanatory environmental variables for both the edaphic and climate models.Model choice was based on Akaike's Information Criterion as a measure of model parsimony (Sakamoto et al., 1986).Pearson correlation coefficients and variance inflation factors were assessed to avoid unnecessary redundancies in the final suite of explanatory variables using the vif.cca function.Models and individual axes were tested for significance using the anova.ccafunction.Vectors were superimposed onto RDA biplots with significance established through 999 permutations using the envfit function, also in the 'vegan' statistical package (Oksanen et al., 2007).We log 10 transformed trait variables of stomatal conductance and rosette size to meet assumptions of homogeneity of variance, and all variables were centred and standardized prior to RDA to provide a comparable unitless scale to compare the trait dataset with climatic and edaphic variables.
To partition the effects of climatic and edaphic factors on trait compositions, we ran three different RDAs: a full model including both climatic and edaphic explanatory variables identified in the forward selection procedure, a partial model of climatic variables while controlling for edaphic effects and a partial model of edaphic variables while controlling for climatic effects.The unique and joint effects of climate and soil for explaining the variation in each morphological and physiological traits were then quantified using variation partitioning (Legendre and Legendre, 2012).The shared variance accounted for by both climatic and edaphic factors was calculated as the difference between the sum of the adjusted R 2 for soil and climate and the adjusted R 2 of the full model encompassing all factors.Marginal significance of individual explanatory variables was established through 999 permutations using the anova.ccafunction in the 'vegan' statistical package (Oksanen et al., 2007).Partial RDAs were used to determine the proportion of variance that was individually constrained by the explanatory variables and the proportion that was 'conditioned' by all remaining variables from which we calculated the proportion of shared variance for each explanatory variable.
The first two components of the PCA describing variation in functional traits accounted for 46.1% of variation across the latitudinal gradient (Fig. 4a).The first axis, which was characterized primarily by a gas exchange gradient, accounted for 26.1% of this variation.The leaf traits that displayed the strongest positive correlations to the first dimension were stomatal conductance (R 2 = 0.91, P < 0.001) and transpiration (R 2 = 0.85, P < 0.001), respectively, while WUE (A/g s ; R 2 = 0.46, P < 0.001) presented the strongest negative correlation to the first axis.The second axis accounted for 20.0% of the trait variation.However, this axis was characterized by
The PCA describing climate cumulatively described 69.2% of variation with the first two principal component axes (Fig. 4b).The first axis, representing local microclimate, accounted for 50.6% of this variation while the second accounted for 18.6%, representing long-term climate averages.The first axis was positively associated with the mean ambient temperature as constrained by the survey period (R 2 = 0.96, P < 0.001), followed closely by local maximum (R 2 = 0.92, P < 0.001), and negatively associated with soil moisture variability (R 2 = 0.92, P = 0.001) and rainfall events as constrained by the survey period (R 2 = 0.91, P = 0.001).The second axis was positively associated with MAT (R 2 = 0.83, P < 0.001) and negatively associated with MAP (R 2 = 0.81, P < 0.001).
The first two components of the PCA describing variation in the edaphic environment across the gradient accounted for a cumulative 72.3% of variation (Fig. 4c).The first axis component, which reflects a gradient of macro-and micronutrient availability, accounted for 53.5% of this variation.Sulphur (R 2 = 0.94, P < 0.001) presented the strongest pos- itive association with the first component, closely followed by boron, ammonium nitrate, potassium, sodium and phosphorus, respectively (R 2 range = 0.82-0.90,P < 0.001).The second principal component axis accounted for 18.8% of the variation.However, this axis corresponded most strongly to a trace metal and acidity gradient, which was significantly positively associated with soil pH (R 2 = 0.91, P < 0.001) and negatively associated with aluminium (R 2 = 0.88, P < 0.001) and iron (R 2 = 0.71, P < 0.001) concentrations.

Soil factors
In the RDA for leaf traits and soil characteristics, only two factors, soil pH and phosphorus content, remained in the model after the forward selection procedure.The RDA model constructed for the edaphic drivers explained 9.3% (R 2 = 0.093, adjusted R 2 = 0.080) of the weighted variance (inertia) in the leaf trait data, as constrained by phosphorus and soil pH (RDA1: Pseudo-F 1,139 = 7.89, P = 0.001; RDA2: Pseudo-F 1,139 = 6.44,P = 0.001; Fig. 5b, Table 1a).Analysis of individual terms indicated that both phosphorus (adjusted R 2 = 0.042, F 1,140 = 8.69, P = 0.001; Table 1b) and soil pH (adjusted R 2 = 0.041, F 1,140 = 8.49, P = 0.001; Table 1b) accounted for a similar amount of variation in the leaf trait data.

Variance partitioning
Variance partitioning showed that overall, climate factors, as determined by the forward selection procedure, explained more variation (adjusted R 2 = 0.192; Table 1a) in the collective matrix of leaf traits than did soil (adjusted R 2 = 0.080; The proportion of variance explained by each axis appears in parentheses, followed by the proportion of inertia explained.R 2 and adjusted R 2 ( adj R 2 ) are overlayed in grey text.The climatic and edaphic factors used in these analyses were the ones revealed to be the most relevant following a forward selection procedure using the ordistep function in the R statistical package 'vegan'.Vector labels correspond to mean atmospheric VPD (kPa), mean maximum temperatures ( • C, Ta max ), mean ambient temperature ( • C, Ta mean ), MAP (mm), RH (%), soil pH (pH) and phosphorus content.

Table 1a
).A partial model of climatic variables while controlling for edaphic effects captured more variation (adjusted R 2 = 0.134; Table 1a) in the plant trait dataset than did the partial model of edaphic variables while controlling for climatic effects (adjusted R 2 = 0.022; Table 1a).Climatic and edaphic factors had a shared explained variance of 5.8% (adjusted R 2 = 0.058; Fig. 6a, Table 1a).Partitioning out the effects of climate, soils and their joint effect for each independent plant trait showcased between 15.1% and 97.6% of the explained variation observed within each trait was captured by the independent effects of climate, and up to 38.5% of the variance was captured by the independent effects of soil (Fig. 6b).Joint contributions ranged between 2.02% and 84.9% (Fig. 6b

Table 1: a)
Results of RDA on trait compositions of S. hispidum Lindl.Models partitioned inertia (variance) into that encompassed by climate and soil factors, as well as a partial model of climatic variables while controlling for edaphic effects (Climatic | Edaphic), a partial model of edaphic variables while controlling for climatic effects (Edaphic | Climatic) and a joint model capturing the shared inertia between both climatic and edaphic models (Climatic ∩ Edaphic) with asterisks denoting model significance ( * , P < 0.05; * * , P < 0.01; * * * , P < 0.001, NA = Not Available).b) Partitioned inertia by simple RDAs of independent climatic and edaphic factors.Proportion constrained corresponds to the variance constrained by the sub-models (in a) and independent factors (in b) relative to the proportion constrained by the full model.Asterisks denote marginal (Type III) effects of individual explanatory variables in a model with all other terms accounted for.

Discussion
Examining where morphological and physiological trait diversity exists, and the environmental factors associated with ITV, can provide insight into how plants are performing in current conditions and their resilience to environmental change.Our research details the physiological and morphological variation within and between five spatially distinct S. hispidum populations and the associations with local climatic and edaphic conditions.There were marked differences in the thermal and hygric conditions along a latitudinal gradient, confirming that northern sites were hotter and drier than southern sites.We found that northern populations presented significantly reduced gas exchange, rosette size and floral investment compared to the southern populations, despite the similar patterns of trait variability (QCD) across populations.
Trait expression in the northern populations is likely a result of morphological and physiological strategies that work to cope with heat and drought stress.Overall, climate factors, and in particular, temperature, accounted for the greatest proportion of variation in independent morphophysiological traits, as well as multivariate trait composition.

Patterns of trait variation
Conforming to our hypothesis, that trait means would vary along the gradient, we found northern populations to present trait means reflective of water conservation strategies such as reduced stomatal conductance, smaller rosettes and limited floral investment compared to the southern populations.Generally, these patterns of form and function across the latitudinal gradient were captured by two main axes of variation (Fig. 4a), with the first axis broadly reflecting a gasexchange gradient and the second reflecting a productivity and morphological gradient.Such axes have emerged in various forms across several previous studies.For example, both Díaz et al. (2016) and Joswig et al. (2022)   ( Joswig et al., 2022).Therefore, morphophysiological trait data may be used to generate predictions of how species may respond to changing environments to gain insights into the mechanisms driving trait-environment relationships.
Intraspecific variation is a key component of functional trait diversity, and adaptive plasticity is a major mechanism by which plants respond to abiotic and biotic stress (Gratani, 2014).High ITV is considered to promote species acclimation, increase competitive capacity and population persistence in dynamic environments and may be a key index for predicting population and species resilience in the face of a changing climate or increasingly degraded landscapes (Valladares et al., 2007;Valladares et al., 2014).However, we found that trait variability was evenly distributed along the latitudinal gradient with all populations presenting statistically indistinguishable differences in the QCD.Consequently, our results do not support the prediction that northern populations would reflect reduced trait variability, although this finding is consistent with previous studies that report no evidence to suggest that there is population differentiation in intraspecific trait variability (Matesanz et al., 2012;Kühn et al., 2021).While we do not distinguish between the trait variation that arises from heritable genetic differences or by phenotypic plasticity, the similarities in trait variability among populations may suggest that there is no differential selection for plasticity across the latitudinal gradient, or there is insufficient heritable variation to drive population differentiation in trait variability for S. hispidum (Goldstein and Ehrenreich, 2021).It could also be suggested that the similar levels of population level variability are the result of similar degrees of heterogeneity, despite different means, in abiotic conditions at different sites (De Kort et al., 2020).However, these conclusions are context dependent and depend largely on the species, environmental heterogeneity and the extent of a species distribution.

Drivers of latitudinal variation
Generally, climatic and edaphic gradients are associated with significant differences in functional traits, relating to the morphology, physiology and reproduction of plants, and they are recognized as important predictors of population persistence in the face of environmental change (Violle et al., 2007).We found marked differences in the thermal and hygric conditions along the latitudinal gradient, such that northern sites were 1.6 • C warmer than southern sites, on average, with maximum temperatures reaching 3.3 • C higher, while there was a 99% increase in mean rainfall from the driest northern site to the wettest southern site during the survey period alone.Moreover, we observed significantly higher macronutrient availability in the second most northern site, John Forrest.While the effects of climate on edaphic factors were not examined here, the patterns observed in soil factors loosely correspond to Albrecht's conceptual model (Huston, 2012), such that soil pH, exchangeable cations and nutrient availability track a unimodal pattern from being low at the extremes for dry and wet conditions (Huston, 2012;Maire et al., 2015).
RDA, together with variance partitioning, allowed the determination that climate factors consistently accounted for a significantly greater portion of the weighted variance in plant traits than that of edaphic factors.Of the climate factors defined within this study, only ambient temperature, maximum temperature, atmospheric VPD, moisture availability and rainfall presented significant associations with variance in plant functional traits along the latitudinal gradient, while only soil pH and phosphorus emerged as significant predictors out of the available soil factors.Separately, the climate factors alone constrained 22.1% (0.192 adjusted R 2 ) of the total inertia in plant traits, while the soil factors constrained 9.3% (0.080 adjusted R 2 ).This result contrasts with findings recorded at global scales, whereby soil predictors were found to explain more variance in individual leaf traits than climate (Ordoñez et al., 2009;Maire et al., 2015).However, this finding highlights the dominant role that climate can play in structuring trait expression in S. hispidum and suggests that future climate change scenarios may lead to trait shifts across populations in response to the effects of warming and drying.Moreover, there is uncertainty surrounding the extent of abiotic stress that this species can tolerate, highlighting the importance of characterizing physiological thresholds and resilience at biologically meaningful scales.
Despite edaphic factors being less correlated with plant traits than climate in our model species, we found that soil pH and phosphorus content explained 4.1% and 4.2% of total inertia, respectively.These results conform, in part, to those observed in a global, interspecific analysis conducted by Maire et al. (2015), finding the same two soil variables to have the strongest explanatory power over the variation in plant traits.While we found the explanatory power of soil pH and phosphorus availability to be substantially less than climatic variables, both soil pH and phosphorus content play key roles in shaping plant growth, resource acquisition and performance (Malhotra et al., 2018;Neina, 2019).When considered across a broader gradient of soil ecotypes, increased alkalinity will often be associated with a higher availability of nutrients and reduce the acquisition costs of key macronutrients, such as nitrogen (Maire et al., 2015).Increased acidity, however, can compromise root form and function and, as a consequence, inhibit CO 2 assimilation (Long et al., 2017), reduce water uptake (Bian et al., 2013) and induce oxidative stress (Martins et al., 2013).We observed a unimodal distribution of phosphorus content along the latitudinal gradient with intermediate populations presenting significantly greater phosphorus availability than that of marginal populations.Photosynthesis is generally reduced in phosphorus deficient soils, although there is some degree of adaptation to be expected for plants growing in environments with reduced phosphorus availability (Lambers et al., 2015).Nevertheless, the unimodal patterns in photosynthesis do correspond to those of phosphorus availability, likely contributing to spatial variance observed in photosynthesis in this study.
The tendency for plants to exhibit morphological and physiological shifts along climate gradients is well known and in accordance with theory to predict shifts in gas exchange and water relations as strategies to minimize water loss and maximize water uptake whilst maintaining physiological functioning (Reich et al., 2003;Westoby and Wright, 2006).We found that the strongest climate driver of variation in the leaf trait data was temperature, with mean ambient and maximum temperatures independently explaining 11.8% and 12.6% of inertia, respectively.This conforms to global scale trends in plants traits as they relate to temperature and precipitation, with Moles et al. (2014) finding that temperature is significantly correlated with 15 out of 21 plant traits while precipitation correlated to a lesser extent, with only six of these traits.Above-optimal temperatures and drying environments can impair plant growth, development, form and function (Bykova et al., 2012), often exerting greater influences in combination than either abiotic stress does independently (Zandalinas et al., 2018).Collectively, the patterns observed in plant traits observed in the northern populations can be seen as strategies to manage the economics of gas exchange and be under stronger selection pressures from climatic factors to mitigate or cope with warming and drying conditions.However, we are limited in our capacity to conclude whether the patterns of trait means are truly reflective of genetic differentiation supporting local adaptation.Therefore, future studies should seek to disentangle the relative contribution of environmental, maternal and genetic effects on plant responses (e.g.growth rates and physiological function) to combined heat and drought stress under controlled "common-garden" scenarios.

Implications for conservation
While our results show short-term responses to interacting environmental drivers, there may be long-term implications for local extinction if trait modifications that impact plant fitness fail to keep up with rapidly changing environments.In the southwest region of Western Australia, daily maximum temperatures are expected to increase up to 2.4 • C by midcentury and 4.2 • C by the end of the century (Jakob et al., 2012;CSIRO and BOM, 2015).Similarly, MAP is expected to decline by 17% by mid-century and 46% by the end of the century under an intermediate-emission scenario (Jakob et al., 2012;Sudmeyer et al., 2016).In this regard, it is possible that by the end of the century, the southern populations may present a phenotype more reflective of those currently observed in northern populations, as the latitudinal gradient represents a thermal gradient difference of 1.4 • C in the warmest quarter, and 39.4% decrease in MAP.The consequences for ongoing warming and drying in the northern populations are less clear, though they may be detrimental if the adaptive capacity of the species is exceeded by the rate of change, or if tolerance thresholds for abiotic stress are surpassed.In this instance, hotter and drier conditions may result in range contractions, seeing the local extinction of northern-most populations, or potentially range shifts or expansions further south to habitats that foster milder conditions (Fitzpatrick et al., 2008), especially considering that climate-related local extinctions are already widespread across the globe (Wiens, 2016).However, to date, there are no studies to report on the intergenerational effects of abiotic stress on Stylidium, nor any specific parameterisation of drought or thermal tolerance.Therefore, it is difficult to predict with greater certainty as to how this species will perform in future climates, highlighting a major gap in our scientific knowledge for not only the focal species of this study, but also the exceptionally diverse flora of southwest Australia.
Unless appropriate conservation and management initiatives are implemented to address the potential threat of climate change, the future remains unclear for northern marginal populations of S. hispidum in warming and drying southwest Western Australia.From a conservation perspective, ex situ storage of seed material may temporarily safeguard population genotypes and serve as a repository source for in situ conservation, restoration or translocation efforts (Merritt and Dixon, 2011).Additionally, it is important to monitor changes in population abundance and trait compositions as climate-associated changes unfold.If vulnerabilities are exacerbated in the future, targeted gene flow (Hufford et al., 2012;Hufford et al., 2016), managed relocation (Butt et al., 2021;Onley et al., 2021) and microclimatic amelioration (Soliveres et al., 2011;Brigham and Suding, 2023) may be considered to facilitate adaptation or open migration pathways, especially for priority species.Nevertheless, it is important to note that the findings, interpretations and speculations posed in this study do not necessarily translate other species, as the results are largely species specific and context dependent, with other trends potentially emerging at different spatial and temporal scales than those observed here.As such, ongoing conservation of at-risk species calls for investment into evidence-based insight into the intraspecific variation in tolerance to interactive stressors associated with environmental change.

Conclusions
Here, we have demonstrated that while intraspecific trait variability may be consistent among populations spanning a latitudinal gradient, patterns of trait means are strongly associated with climate, and to a lesser extent, edaphic factors.While the trait-environment correlations identified in our study should not be assumed to be direct causal relationships, the dominant impact of climatic factors shaping functional trait expression along the gradient suggests that trait shifts will occur with future climate change.Disentangling the abiotic drivers of functional trait variation will be key to understanding the mechanisms underpinning local adaptation and population-level responses to ongoing environmental change.Understanding these dynamics is of global relevance; however, it is particularly important within the context of our study system, where increased temperatures and declining rainfall risk the 'transition or collapse' of the vulnerable Jarrah Forest ecosystems (Lawrence et al., 2022).

Figure 1 :
Figure 1: Conceptual diagram showcasing the hypotheses for change in intraspecific trait variation (ITV), calculated as the quartile coefficient of dispersion (QCD), and functional trait means for S.hispidum populations along the local latitudinal gradient.We expected greater ITV (represented by green arrows) in southern populations and for trait means to vary along the gradient (represented by green dots), with northern populations having trait means reflective of water conservation strategies (e.g.reduced stomatal conductance, smaller plant sizes and fewer flowers).Please note that not all functional traits were expected to exhibit the visualized positive trend across the latitudinal gradient (i.e.WUE was expected to decrease from north to south); the linear trend line is merely indicative of differences in means, irrespective of direction and slope.

Figure 2 :
Figure 2: Location of study sites (♦) and climatologies distributed across the ∼160-km latitudinal gradient along the Darling Scarp in the southwest corner of Western Australia (inset) with current distribution of S. hispidum populations (•) overlayed on the Interim Biogeographic Regionalisation for Australia with the NJF bioregion shaded in green.Monthly climate profiles composed of mean monthly precipitation (blue bars) and mean maximum monthly temperature (black dashed line) are constructed from the nearest available weather station data accessed from the Bureau of Meteorology, Perth, Western Australia.MAP (mm; blue text), mean temperature in the warmest quarter ( • C; black text) and elevation (m; green text) are inset within each graph with precipitation and temperature data extracted from BioClim.
Figure 3: a) Spatiotemporal variation in functional traits across at least 10 S. hispidum plants (mean ± standard error) for each site along the surveyed latitudinal gradient by sampling month.Error bars, from left to right, correspond to sites from north to south-Julimar State Forest, John Forrest National Park, Bungendore Park, State Forest (Dwellingup), State Forest (Hoffman).b) Intraspecific trait variation as measured by the QCD across the latitudinal gradient.Error bars (mean QCD ± standard error), from left to right, correspond to sites from north (N) to south (S)-Julimar State Forest, John Forrest National Park, Bungendore Park, State Forest (Dwellingup), State Forest (Hoffman).

Figure 5 :
Figure 5: RDA biplots showing plant functional trait composition as constrained by a) climatic factors and b) edaphic factors.Opaque points (±95% confidence intervals) in panels a) and b) represent site centroids, showcasing the latitudinal variation of survey sites, from north (N) to south (S)-Julimar State Forest, John Forrest National Park, Bungendore Park, State Forest (Dwellingup), State Forest (Hoffman).The proportion of variance explained by each axis appears in parentheses, followed by the proportion of inertia explained.R 2 and adjusted R 2 ( adj R 2 ) are overlayed in grey text.The climatic and edaphic factors used in these analyses were the ones revealed to be the most relevant following a forward selection procedure using the ordistep function in the R statistical package 'vegan'.Vector labels correspond to mean atmospheric VPD (kPa), mean maximum temperatures ( • C, Ta max ), mean ambient temperature ( • C, Ta mean ), MAP (mm), RH (%), soil pH (pH) and phosphorus content.